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CN , Abstract 

(— I I A hyperscaling relation for the critical exponents of absorbing phase transitions is tested in the 

O ■ 

D ' bosonic pair contact process with diffusion. To this end spreading is considered, i.e. the time 



evolution out of an initial seed. It is shown that like in the case of spatial homogeneous initial 
conditions the autocorrelation function exhibits a phase transition at the critical point of the first 
moment. Some of the critical exponents can be determined exactly which is an unusual property 
of an absorbing phase transition and provides a possibility to test the hyperscaling relation. In the 



r^ . case of the bosonic pair contact process with diffusion three sets of exponents can be considered 

o, 

O I referring to the number of particles, number of pairs and number of active sites. It is argued that 



in special cases it is generally impossible to produce adequate data numerically. 
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I. INTRODUCTION 

A possibility to investigate non-equilibrium critical phenomena is given by absorbing 
phase transitions. These are transitions from an active fluctuating phase with a finite particle 
density to an absorbing state where any dynamics is suppressed. A central question in the 
theory of critical phenomena is the determination of the universality class of a given system. 
In the field of absorbing phase transitions rather robust universality classes have been found, 
e.g. the class of directed percolation (DP) and the parity conserving universality class (PC). 
For a review see Ref. [l|. 

In order to determine the respective universality class the calculation of critical exponents 
plays a central role. One possibility to determine critical exponents is to consider the scaling 
of stationary quantities with respect to the distance from the critical point A = p — pc where 
p is the control parameter and Pc the critical point. One defines the exponents (3,(3', z/_l, u^ 
by 

(1) 

X± oc |Ars 

XII oc I Aril, 

where p^**^* is the stationary particle density. Poo is the ultimate survival probability, i.e. the 
probability that a randomly chosen site belongs to an infinite cluster, and x±5 X\\ ^ire the 
spatial and temporal correlation lengths. 

Another possibility is given by dynamical scaling where the time dependence of the 
quantities when started from an initial seed is used to define the exponents. At the critical 
point, A = 0, one defines: 

P{t) oc r\ 
p{t) oc t-, (2) 

{N{t)) oc t^ 

where P{t) is the probability that a system survives at time t, i.e. that there are still active 
particles left, p{t) is the particle density inside an active cluster and {N{t)) is the particle 
number averaged over all, i.e., active and inactive, systems. This set of exponents is not 



independent of the previous set, one can deduce generally for second order transitions |l| 

<5 = /37z/||, a = /3/i/||. (3) 

Furthermore the following argument gives another relation between the exponents: The 
particle density inside an active cluster is given by the particle number in a specific active 
system, N{t) divided by the spreading region R{t) 

Pit) = (^)active, (4) 

here (•)active indicates that the average is taken only over active systems. For large times 
one expects that 

as the spreading region scales like t^^^ where d is the dimension and z = i^\\/i^± is the 
dynamical exponent. Thus we get the hyperscaling relation 

0-d/z = -{a + S). (6) 

This rather intuitive derivation should hold for first 2] and second order phase transitions. 
A more detailed derivation of the hyperscaling relation for second order transitions can be 
found in Refs. [l|, y]. It has also been shown that a hyperscaling relation can be defined for 
the case of coupled systems, see Ref. j^ and references therein. 

As there is hardly any model used for investigating absorbing phase transitions that can 
be solved analytically, all critical exponents have to be determined numerically. An exception 
to this is the bosonic pair contact process with diffusion (bosonic PCPD) for which a field 
theoretic approach due to Howard and Tauber [5] is available. In this context the term 
'bosonic' refers to the property that there is no exclusion rule for the particles, each lattice 
site may be occupied by any number of particles. This leads to a theoretical description in 
terms of bosonic operators instead of fermionic operators in the case of exclusion models. 

A drawback of this model is that it exhibits a first order transition and is thus not suitable 
for deciding the universalitv class of the PCPD with particle number restriction, which is 
still an open problernja, [O, Isl, Isl lll| • For a comprehensive review of the current state of the 
art we refer to Ref. jll|. 

However, an exceptional property of the bosonic PCPD is the fact that some of the 
critical exponents are known exactly ^]. This information provides a possibility to test 



the commonly considered hyperscaling relation Eq. (jH)) for the critical exponents. As still 
some quantities are not accessible analytically Monte Carlo simulations are used to get the 
complete set of exponents. It turns out that in some cases simulations of the bosonic PCPD 
are misleading in general. 

II. MODEL AND FORMALISM 

We consider the following process as introduced in Ref. p|: On a infinite d-dimensional 
cubic lattice particles ('^4') are diffusing with rate D, in each spatial direction. Additionally 
they branch and annihilate: k > 1 particles A are created with rate fi out of any set of 2 
particles, and / G {1,2} particles are annihilated with rate A out of any set of 2 particles: 

2A A {m + k)A 



2A ^ {p-l)A 



A- ^ -A. 



(7) 



The number of particles on each lattice site is not restricted - the creation and annihilation 
processes take place on one lattice site. Thus the bosonic representation of the process is 
used. One special case is the PCPD, where / = 2 and k = 1. 

For Xl > fik the particles die out according to a power law, while for XI < fik the particle 
density diverges. In analogy to the exclusion model we call the rate which divides the two 
different behaviors the "critical" rate, 



Ac = fj^k/l 



for given n 12] . For this rate the total particle number 
For the details of the formalism we refer to Refs. |5 



(8) 

is constant for all times. 
12| . The time evolution of the particle 



annihilation operators a(x) is derived for the case A = /i/c/Z: 

|(a(x)) =D J2 {(«(x - k)) + (a(x + k)) - 2(a(x))} 
fc=i 

-(a(x)a(y)) = /^ J^ { (a(x)a(y - k)) + (a(x)a(y + k)) + 

A;=l 

(a(x - k)a(y)) + (a(x + k)a(y)) - 4(a(x)a(y)) } 



(9) 



|((a(x))^) =2D J2 { («(x)a(x - k)) + {a{^)a{^ + k)) - 2(a(x)2)} 
fc=i 

+ /iA;(A; + Z)(a(x)'") 

where k = k(/c) = (. . . , 0, 1, 0, . . .)^ is the A;-th unit space vector. The particle density is 
given by p(x, t) = (^(x)) = (a(x)) and the correlation function by (n(x)n(y)) = (a(x)a(y)) + 
5x,y(a(x)). 

III. ACTIVITY SPREADING 

In the theory of absorbing phase transitions beneath spatial homogeneous initial condi- 
tions often the following scenario is used: Initially the lattice is empty except for the origin 
where just as many particle are located as needed for the dynamics to start. It is then inves- 
tigated how this activity spreads into the system. In this section the question is addressed 
what can be learned from this initial condition in the bosonic PCPD. 

After presenting analytically solvable cases we test the hyperscaling relation of the critical 
exponents. To this end additional information is needed which can be obtained from Monte 
Carlo simulations. 

A. Analytical calculation 

As initial condition of the system we choose two particles at the origin 

p(x,t = 0)=2(5x,o. (10) 

For the average particle number (X]x'^(x)) we recover the same result as for the density 
in the case of spatially homogeneous initial conditions: for XI < fik the particle number 
diverges ("active phase") while all particles die out for \l > ixk ("inactive phase"). In 



the active phase a spreading cone forms, i.e. a growing region with non-zero density. If 
A = ^k/l is chosen the time evolution of the density is simply a lattice diffusion equation 
and the solution is given by 

p(x,t) = 2e-2'^%,(2t)-...-4,(2t). (11) 

For large arguments x and t this function asymptotically approaches a Gaussian distribution, 
thus the dynamical exponent is z = 2 and the average particle number is constant, 6* = 0. 
By defining 

Fx(r,t) = (a(x)a(x + r)) (12) 

and rescaling time by 



the time evolution of the second moment can be rewritten as 

d 

2 



1 '^ 
F^{r, t) =- J2 {^x(r - k, t) + F,(r + k, t) 



k=l 



+F,_k(r + k, t) + F,+k(r - k, t) - 4F,(r - k, t)} 

+ (5r,oFx(0,t), 



(14) 



where 



21) ^ ^ 

We rescaled time by the factor 1/(2D) instead of 1/D in order to keep the notation consistent 
with the previous publication, Ref. |l2l |. 

Equation (fT^ can be solved using a two-component Fourier-transformation 



/(s,q,t)=5^5^e--e-^'i^F.(r,t), 



X r 



(16) 



r d sd a 

Fx(r,t)=y ^^e-V-/(s,q,t). 

The differential equation for the Fourier-transform / can be cast into the form 

-/(s, q, t) = -h{s, q)/(s, q, t) + «F(s, t), (17) 

with the dispersion relation t>(s, q) = — J2k=i (cos(gfc) + cos(gA: — Sk) — 2) and 

F(s,t) = 5^e— F.(0,t) = J ^/(s,q,^)- (18) 



Integration yields 

/(s, q, t) = e-^^^''^)* |/(s, q, 0) + a £ drF(s, r)e''(^''')^ | . (19) 

The initial condition equation pOj) reads in Fourier space /(s, q, 0) = 2. Thus we get for 
the correlation function 

Fx(r,t) = 2A(x,r,t) + aV f drF,,(0, r)A(x - x', r, t - r), (20) 

where 

./ X f d'^sd'^q / / X . . X ,„.x 

-4(x,r,t)= / — — ^exp(-f(s,q)t + zsx + 2qr). (21) 

An analytical solution of this integral equation could not be found but for the sum of the 
auto correlations , 

F(0,t) = 5^F.(0,t), (22) 

X 

the situation simplifies because of the following identity 

bit) := J]A(x-x',0,t) = J]A(x,0,t) = J (|^e-(°''')* = (e-^*J(2t))'^ (23) 

where / is a modified Bessel function. We thus get 

/•* . 
F{0,t) = 2b{t)+a drF(0,r)fe(t-r). (24) 

Jo 

This equation is exactly the one for the Lagrangian multiplier in the mean spherical model 

and we can use the already known results here [l3|. Using a Laplace transformation, 



we get 



F{p) = / dte~P*F(0,t), (25) 

Jo 



m = -^. (26) 

1 — ao[p) 



An analysis of the behavior for small p gives the late time behavior of -F(0, t). Depending 
on the dimension we find a phase transition with respect to a. The critical point is given 
by the same Uc as found before for spatially homogeneous initial conditions as calculated 
in Ref. [I4I. Above the critical point, a > ac, the sum of the autocorrelations diverges as 
before. At the critical point a = ac the sum of the autocorrelations follows a power law 
i-(2-'i/2) for 2 < d < 4 and approaches a constant for d > 4. Below the critical point a < ac 



the sum of the autocorrelations follows a power law t~'^/^. Not that for c? = 1 the critical 
point is zero, ac = |l2l |. 

The behavior below the critical point can be understood using the interpretation from 
spatial homogeneous initial conditions that in this regime the diffusion is dominant. A 
diffusive system without reactions shows the same late time behavior of the autocorrelator: 

((n(x,t))^) = ^raV(x,t) 

n 

= 5^n2(P(x,t|O,0)r 

n 

X 

= EE^'K''*^^.(2^)-----^^.(2^))" (27) 

X n=2 



$^4y" -|^exp (z(q + q')x - (w;(q) + ^(q')) t) + 



t— >oo 



4 (Svrt) 



-d/2 



where Pn(x, t) is the probability to find n particles at site x at time t which for independent 
particles starting from the origin can be expressed as products of the propagator P(x, t|0, 0). 
Above the critical point the reaction processes are dominant. As on average the particle 
number in each system is 2 in most of the systems the particles have to vanish in order that 
in few systems a divergence of the second moment is possible. This has crucial influence on 
the possibihty to simulate the process as discussed in the next section. 

B. Hyperscaling relation 

1. Analytical predictions 

The hyperscaling relation Eq. (jH)) shall now be verifled in the bosonic PCPD where we 
know after all some of the exponents exactly. The arguments given for the hyperscaling 
relation should hold irrespectively of the type of density/number which is measured. While 
in the description of the process with exclusion interaction there are only two possibilities 

8 
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FIG. 1: Time evolution of the average number of active sites for different parameters, the reaction 
rates are fixed to /i = 2 and A = 1. For d = 1 always a > ac, for d = 3 and d = 5 the diffusion 
constant D is chosen such that a < ac, a = ac and a > ac (with decreasing D). One observes 
that only in one dimension the number of active sites does not approach a constant, a fit yields 
^sites ^ _o.450. 
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FIG. 2: Time evolution of the average number of particles, parameters as before. As expected 
from the analytical calculations the number of particles is constant for any set of parameters. 



- the number of particles and pairs, i.e. two particles on neighboring sites - in the bosonic 
description there are three possibilities. The number of active sites, these are sites with 
at least one particle, have to be distinguished from the number of particles. So one may 
consider the number of active sites, of particles or of particle pairs. The number of particles 
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FIG. 3: Time evolution of the average number of pairs, parameters as before. In contrast to the 
analytical result one can not observe the divergence for a > ac as explained in the text. For 
a = ttc and a < Oc the fluctuations are still very large and the agreement of the slopes with 
the analytically predicted values is not very high. Fitting yields ^^^3*^ = —0.23 to be compared 

—1.7 (analytically: —1.5), 



to -0.5, ^n' 



-0.76 (analytically: 0) for a = Oc and ^^^3^ 



J^^s" ^ -2.5(±0.5) (analytically -2.5) for a < Uc- 



is constant at the critical point, thus we have 6'p'*''* = 0. The number of pairs is given by 
^^((a"''(x))^(a(x))^) = X]x('^^(^)) whose behavior is calculated in subsection IIII Al 



] pairs 



1 



—d/2 a < ac 
-(2 - d/2) a = ac and 2 < d < A 
a = ar and d > 4. 



(28) 



For a > ac the number of pairs diverges exponentially and thus no exponent can be defined. 

In (i > 2 it is expected that 6 = 0: After some initial time there is a fraction of systems 
consisting only of two particles at different sites. They are diffusing freely on the lattice and 
hence their distance vector describes as well a random walk in d dimensions. As a random 
walk is transient for d > 2 the probability that the two particles ever meet again is zero and 
consequently those active systems survive for ever and one concludes 6 = 0. 

Additionally we have shown exactly that z = 2. The remaining exponents, 
^sites^ ^sites^ aP'^'^*, a^'^"^ and 6d=i have to be determined numerically in a Monte-Carlo simu- 
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FIG. 4: Time evolution of the survival probability, parameters as before. As predicted, only for 
d = 1 the curve does not approach a constant, a fit yields in this case 5^=1 = 0.471. 

lation. To sum up, the theoretical considerations predict (see Eq. (0) and Tab.Q}: 
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a = ac and 2 < rf < 4 

a = ac and d > 4 



(29) 



2. Simulations 

For the Monte-Carlo (MC) simulation of the bosonic PCPD a list of active sites i is used. 
In contrast to the model with exclusion interaction the number of particles rii on the sites 
has to be tracked. The number of pairs on each site is then given by ni{ni — l)/2. For 
each system in one MC time step DY^-rii diffusion processes, fiYli^ii^i ~ l)/2 creation 
and A^^nj(nj — l)/2 annihilation processes take place on average. One of these possi- 
bilities is chosen randomly according to its statistical weight and the time is updated by 
t ^ t + [-D ^^ rij + (/i + A) ^j Tiiirii — l)/2]~ . A difficulty is that the number of possible 
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1000 



FIG. 5: Time evolution of the density of active sites, parameters as before. By fitting these curves 
we observe that the exponents are given by a^^^^^ = d/2 in good agreement, for exact values see 
table HH 

processes varies extremely from system to system as the number of pairs fluctuates enor- 
mously. Consequently it is not convenient to determine a target time and simulate each 
system up to this time one after each other because for a badly estimated target time the 
program might get stuck in only one of the systems with a large number of pairs. We rather 
determined target times in small steps up to which the systems were simulated step by step. 
Although the systems have to be kept simultaneously in memory one gains the advantage 
that the results can be tracked during the simulations and one does not have to estimate 
the maximal simulation time in advance. 

Still the simulation of the process takes much effort. Therefore compared to the simula- 
tions of the model with exclusion interaction only small times could be simulated. Especially 
in the case of diverging autocorrelations reliable results are computationally demanding. 
Standard simulation methods simply fail in this case as the number of needed systems in 
the ensemble is far too large. This number of systems M can be estimated as follows: The 
average number of pairs diverges as (^N^^"'^) ^ A^Q^''^*'exp(t/r) while the average number of 
particles A'"o is constant. To determine a lower bound for M one may assume that all the 
particles available in the simulation MNq pile up at one site of a single system. Then the 
number of pairs is well approximated by {MNqY /2. This has to be equal to MNq'^"^ expit/r) 
as this is the only contribution to the ensemble average of the number of pairs. Thus we 
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FIG. 6: Time evolution of the density of particles, parameters as before. By fitting we again observe 
that the exponents are given by aP'^'^* = d/2 in good agreement, for exact values see table Hll 



conclude that 



^pairs 



^^0 



(30) 



systems are needed in order to allow for the divergence of autocorrelations. 

In our simulations we have typically r ^ 1, for example for the 3(i-case with D = 
0.5, /i = 2, A = 1 one gets r = 0.86. Already for the simulation time t = 100 one would 
need an ensemble consisting of roughly 10^° systems to observe the divergence numerically. 
Consequently it is not expected that for a > ac (which is especially true in d = 1) the 
simulation produces correct results. 

Simulations with several parameters were performed: As dimensions d = 1, d = 3 and 
d = 5 are chosen in order to simulate systems below the lower critical dimension, above 
the lower but below the upper critical dimension and in the mean-field regime of the phase 
transition of the second moment. The reaction rates are fixed to /i = 2 and A = Ac = 1 
such that the average number of particles remains constant in time. For d = 3 and d = 5 
the diffusion constant D is chosen such that a < ac {D = 10 for d = 3 and D = 1 ioi d=5), 
a = ac {D = 0.757 for d = 3 and D = 0.346 for d = 5) and a > ac {D = 0.5 for d = 3 and 
D = 0.1 ioT d = 5). In the ensemble M = 10^ systems are simulated in parallel. 

Fig. ^ shows the time evolution of the average number of active sites. For d = 1 the 
number of active sites decreases according to a power law with an exponent 6'^=*f = —0.45. 
For d = 3 and d = 5 for all parameters the curves approach a constant. As explained above 
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this is expected as a certain fraction of the systems will consist of at least two particles 
performing a random walk without meeting again. 

In Fig. 121 the time evolution of the average number of particles is shown. In agreement 
with the analytical calculations the particle number is constant in time and thus 6^^^^ = 
for all parameters. 

The time evolution of the average number of pairs is shown in Fig. 121 Among the three 
quantities considered this is clearly the most fluctuating one. The simulation indeed fails 
in reproducing the exponential divergence of the number of pairs for a > Oc as discussed 
above. While for d = 1 the number of pairs is increasing according to a power law, for the 
higher dimensions it approaches a constant. Interpreting the results for a > ac has thus to 
be done carefully. 

The survival probability P{t) is shown in Fig. El It is verified for d = 3 and d = 5 that 
P{t) approaches a constant while for d = 1 it decays according to a power law. The fitted 
values for 6 can be found in table IHl For d = 3 and d = 5 the values are very close to zero 
and 6d=i ~ 1/2. 

The directly measured densities are shown in Fig. El (active sites) , Fig. IHl (particles) and 
Fig. [71 (pairs). For the densities of active sites and particles the fitted values for a are in 
most cases in good agreement with the values obtained by the hyperscaling relation, only 
for d = 1 and d = 5,a > ac larger deviations appear. As most of the exponents take simple 
values, we conjecture in these cases that q;^=\^ = 1/2 and a^=f = 5/2 and that the measured 
deviations result from the numerical problems. An obvious disagreement between measured 
values and the hyperscaling relation is found for a"^^^ , where the directly measured value is 
approximately 1/2 while the hyperscaling relation predicts it to be approximately zero. The 
question arises whether the hyperscaling relation is violated or which of the values is wrong. 
As the hyperscaling relation turns out to hold in the other cases we believe it to hold in this 
case as well - the inacuracy in the MC results stems from the fact for d = 1 the numerical 
problem is always present because Oc = 0. We conjecture that 6 = 1/2 and a'^^l = 0. This 
can be imagined as follows: In the surviving systems the active regions spreads diffusively 
and inside the active region the reaction kinetics generates a constant density. Consequently 
the particle number increases in these systems proportional to t°'^, but as more and more 
systems die out according to t~^'^ the particle number averaged over all systems remains 
constant. In rf = 3 and d = 5 the situation is different, the particle density inside the active 
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region decays and the number of active systems remains constant. 

The values for a^'^^'^^ for a < ac could not be determined with high accuracy due to 
the high fluctuations. It would be surprising if the hyperscaling relation did not hold for 
this quantity but the accuracy of our data allows neither for proving nor for disproving the 
relation in this case. 

IV. CONCLUSIONS 

The evolution of the system out of an initial seed (activity spreading) is investigated 
in the bosonic pair contact process with diffusion. It is shown that the second moment 
exhibits a phase transition with the same critical point as for spatially homogeneous initial 
conditions. Above the critical point the sum of autocorrelations diverges and below the 
critical point it decreases according to a power law. This power law behavior can be related 
to purely diffusive dynamics. This shows that below the critical point it can be neglected 
that the particles react because most of the time the particles are diffusing freely. The time 
during which two particles occupy the same lattice site is too short to react because below 
the critical point diffusion dominates above reactions. 

We tested a hyperscaling relation for the dynamical critical exponents. To this end 
exponents have to be determined which are not accessible analytically and are thus calculated 
numerically in a Monte-Carlo simulation. It is shown that for the case of a diverging second 
moment it is impossible to produce accurate data as the necessary size of the ensemble 
diverges exponentially in the desired simulation time. At or below the critical point of the 
divergence of the second moment the hyperscaling relation can be verified. It turns out that 
good accuracy can be achieved for the number and density of active sites and the number 
and density of particles while measuring these quantities for the number of pairs is difficult 
due to large fluctuations. 
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TABLE I: The exactly calculated exponents. Not defined exponents are represented by '-' and 
the exponents to be determined in the simulations with '?'. The exponents 0^^'^^^^ a^^^'''^, a^^''* and 
^pairs j^g^yg ^Q )2g determined numerically. 




1000 



FIG. 7: Time evolution of the density of pairs, parameters as before. Due to the fluctuations 
the accuracy of the exponent a^^^'^^ is not very high. A fit yields the values a"^^^^ = 1.92 to be 
compared to the theoretical value 2, a^^'5^ = 3.76 to be compared to 2.5 for the case a = ac and 
^pairs _ 2 gg (analytically: 3), 0^=5*^ ~ 6 (analytically: 5) for the case a < ac- 



*d=3 
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d 


D 


,5 


/Qsites 


Q,sites 


sites 
"calc 


Q,part 


"calc 


Q,pairs 


pairs 
"calc 


1 


5 


0.471 


-0.450 


0.574 


0.479 


0.517 


0.029 


- 


- 


3 


0.5 


0.009 


0.0108 


1.48 


1.49 


1.49 


1.5 


- 


- 


3 


0.76 


0.006 


0.00429 


1.49 


1.49 


1.49 


1.5 


1.92 


2 


3 


10 


0.0002 


2-10-5 


1.50 


1.5 


1.50 


1.5 


3.58 


3 


5 


0.1 


0.004 


0.15 


2.29 


2.35 


2.28 


2.5 


- 


- 


5 


0.35 


5-10-5 


0.0048 


2.49 


2.5 


2.49 


2.5 


3.76 


2.5 


5 


1 


<io-io 


3-10-5 


2.51 


2.5 


2.51 


2.5 


w6 


5 



TABLE II: The numerically determined values of the exponents compared to the values expected 
form the hyperscaling relation ©. For a > Oc the exponent a'^^^^'^ is not defined. It can be seen 
that for the number of sites and the number of particles the hyperscaling relation is satisfied in 
good agreement while for the number of pairs larger deviations appear. 
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